#include "gtest/gtest.h"
#include "sgEigen.h"

//#include "DataStructure/Common/include/Constant.h"

#include <complex.h>

using namespace SG::Algebra;  // Open the namespace
//using SG::DataStructure::Common::TOLERANCE;

class EigenTest : public ::testing::Test
{
public:
    static void SetUpTestSuite ()
    {
    }
};

/*
* Test the A*x = b given by SuperLU User Guide.
* where A is a 3X3 matrix and b is a 3X1 matrix.
* 
* A =
*   [
*       06.00,  07.00,  05.00;
*       07.00,  13.00,  08.00;
*       05.00,  08.00,  06.00
*   ]
*
* b = [9.0; 10.0; 9.0]
*
* OctaveOnline give the solution as followings,
* x = inv(A)*b = [1; -1; 2]
*
* see https://octave-online.net/
*/
TEST_F (EigenTest, doubleTest)
{
    TEigenLLT<Real_t> solver;

    typedef Eigen::Triplet<Real_t> Triplet;
    std::vector<Triplet> triplets;

    triplets.push_back(Triplet(0, 0, 6.0));
    triplets.push_back(Triplet(0, 1, 7.0));
    triplets.push_back(Triplet(0, 2, 5.0));

    triplets.push_back(Triplet(1, 0, 7.0));
    triplets.push_back(Triplet(1, 1, 13.0));
    triplets.push_back(Triplet(1, 2, 8.0));   

    triplets.push_back(Triplet(2, 0, 5.0));
    triplets.push_back(Triplet(2, 1, 8.0));
    triplets.push_back(Triplet(2, 2, 6.0));

    int m = 3; // row count
    int n = 3; // column count

    SpMat_t mat (m, n);
    mat.setFromTriplets(triplets.begin(), triplets.end());

    solver.compute (mat);

    //Check
    auto info1 = solver.info ();
    ASSERT_TRUE (Info_t::Success == info1);

    // solution from OctaveOnline
    // x = inv(A)*b
    // see https://octave-online.net/
    Real_t exact[] = { 1, -1, 2 };

    Vec_t b (m);
    b[0] = 9.0;
    b[1] = 10.0;
    b[2] = 9.0;

    Vec_t x (m);
    solver.solve (b, x);

    // Check
    auto info2 = solver.info ();
    ASSERT_TRUE (Info_t::Success == info2);

    // for (int i = 0; i < m; ++ i)
    // {
    //     std::cout << exact[i] << ", " << x[i] << std::endl;
    // }

    // for (int i = 0; i < m; ++i)
    // {
    //     ASSERT_NEAR (exact[i], x (i), TOLERANCE);
    // }
}

TEST_F (EigenTest, complexTest)
{
    TEigenLLT<Complex_t> solver;

    typedef Eigen::Triplet<Complex_t> Triplet;
    std::vector<Triplet>                    triplets;

    triplets.push_back (Triplet (0, 0, 6.0));
    triplets.push_back (Triplet (0, 1, 7.0));
    triplets.push_back (Triplet (0, 2, 5.0));

    triplets.push_back (Triplet (1, 0, 7.0));
    triplets.push_back (Triplet (1, 1, 13.0));
    triplets.push_back (Triplet (1, 2, 8.0));

    triplets.push_back (Triplet (2, 0, 5.0));
    triplets.push_back (Triplet (2, 1, 8.0));
    triplets.push_back (Triplet (2, 2, 6.0));

    int m = 3;  // row count
    int n = 3;  // column count

    ComplexSpMat_t mat (m, n);
    mat.setFromTriplets (triplets.begin (), triplets.end ());

    solver.compute (mat);

    // Check
    auto info1 = solver.info ();
    ASSERT_TRUE (Info_t::Success == info1);

    // solution from OctaveOnline
    // x = inv(A)*b
    // see https://octave-online.net/
    Complex_t exact[] = { 1I, -1I, 2I };

    ComplexVec_t b (m);
    b[0] = Complex_t (9.0I);
    b[1] = Complex_t (10.0I);
    b[2] = Complex_t (9.0I);

    ComplexVec_t x (m);
    solver.solve (b, x);

    // Check
    auto info2 = solver.info ();
    ASSERT_TRUE (Info_t::Success == info2);

    // for (int i = 0; i < m; ++ i)
    // {
    //     std::cout << exact[i] << ", " << x[i] << std::endl;
    // }

    // for (int i = 0; i < m; ++i)
    // {
    //     ASSERT_NEAR (exact[i].real(), x (i).real(), TOLERANCE);
    //     ASSERT_NEAR (exact[i].imag(), x (i).imag(), TOLERANCE);
    // }
}